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1.0 Introduction 


This Interim Report summarizes the work performed under Modification 
7 of Contract NAS8-33691 entitled "Study of Tethered Satellite Active 
Attitude Control." Under this modification, SAO adapted existing software 
for the study of tethered subsatelllte rotational dynamics to Its present 
data processing system, developed an analytic solution for a stable con- 
figuration of a tethered subsatelllte, compared the analytic and numerical 
Integrator (computer) solutions for this "test case" In a two-mass tether 
model program (DUMBEL), modified the existing SAO multiple-mass tether 
model (SKYHOOK) to Include subsatellite rotational dynamics, verified this 
modification with the analytic "test case," and demonstrated the use of 
the SKYHOOK rotational dynamics capability with a computer run showing 
the effect of a single off-axis thruster on the behavior of the sub- 
satellite. 

SAO is now In a position to develop subroutines for specific 
attitude control systems and to apply them to the study of the behavior 
of the tethered subsatellite under realistic on-orbit conditions. Such 
studies could also include the effect of all tether "inputs," including 
pendular oscillations, air drag, and electrodynamic Interactions, on the 
dynamic behavior of the tether. 


2.0 A Simplified Model of the Motion About the Center of Mass of a 

Tethered Subsatelllte 

2.1 Introduction 

The dynamics of a tethered subsatelllte about Its center of mass Is 
quite complex In the general case and required sophisticated changes to the 
SKYHOOK program.. SKYHOOK has been rewritten to deal with any tethered sub- 
satellite mission requiring high accuracy attitude prediction and control 
such as astronomical or earth observation missions or, more generally, any 
mission which depends heavily on Isolation of the package from the noise coming 
from the connection to the Shuttle through the tether. For acquiring experience 
and physical insight into the dynamics of this system and for providing an 
analytical solution against which the program changes could be checked, we 
have studied a simple configuration which can be solved analytically. 

2.2 The Simplified Model 

We consider the shuttle in a circular orbit (nominal low orbit at 150 
n.m. = 270 km altitude here, but this is not a very critical constraint) and 
assume that the orbit is polar or equatorial. The orbital plane is therefore 
invariable in an inertial frame. 

We also assume that the subsatellite is axially symmetric with the point 
of attachment of the tether (0) being a point on the symmetric axis of the sub- 
satellite different from G, the center-of-grayity of the subsatellite. We also 
call C the moment of inertia of the subsatelllte with respect to the symmetry 
axis and A, the moment of inertia with respect to an equatorial axis. In 
this model we neglect the dynamics of the tether induced by the rotational 
motion of the subsatellite and we assume also that there is no significant 
in-pl&ne or out-of-plane oscillation of the tether and subsatelllte. 
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W* neglect also the longitudinal oscillation of the tether-subsatelllte system 
Thft motion of the subsatelllte Is therefore equivalent to the motion of a 
axis-symmetric rigid-body about Its center of mass when a force constant 
In magnitude and rotating uniformly about an axis normal to the force and 
constant in orientation (In this case normal to the orbital plane) Is 
applied to a point, P, of Its symmetry axis. Other torques and the effects 
of the precession of the orbital plane may be Introduced into the computer 
model but are difficult to deal with analytically. 

2.3 Notation and Reference Frame 

In Figure 2-1 E Is the center of Earth; S, the Shuttle considered 
as a point; 0, the attachment point of the subsatelllte; a., the unit vector 
from E to S to 0; c, the unit vector normal to the orbital plane; jb e e x 
the unit vector tangent to the orbit; G the center-of-mass of the sub- 
satellite; d, the distance from 0 to G; k, a unit vector in the direction 
from 0 to G; i and two orthogonal unit vectors in the plane normal to the 
symmetry axisO-Gwhich equals kd; N, the nodal line intersection of the plane 
0, 1, j/with the orbital plane; ij > , 4>, e, the Eulerian angles of the frame 0, 
i» i* L» with respect to the frame 0, a., b, c and £, the gravity gradient 
force acting on G which is equal to the tension in the tether. We neglect 
the gravity gradient torque acting on the body which is negligible if the 
dimensions o* the body are small with respect to the length of the tether (&) 

and which it, true for virtually all subsatellites. In good approximation we 
may write: 


E - 3gb)M(£) [l+o(“) + o 


(C-A)n 2 a 


\mgTiT 


] 
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\ EARTH 



ATTACH POINT (0) 


Figure 2-1. Reference Frame for Motion of Subsatellite About Tether Attach Point 
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Here g(0) is the acceleration of gravity at 0, M is the mass of the subsatellite; 
n, the orbital mean motion, o(x) is a coefficient of magnitude x with respect to 
£» and where the quantities denoted by a are Infinitesimally small. As an 
example, If d ■ 0.5m, i ■ 50 km, C-A ■ M(0.5m) 2 , n ■ 10" 3 , g(0) ■ 9 m/sec 2 , 
we have 


and 


a(f) • 10-5 


(2-1) 



* 0.25 x 10" 6 x 6.7 x 10 6 
0.5 • Sx 10 4 • 9 


- lO" 5 


( 2 - 2 ) 


This evaluation gives an idea of the magnitude of the terms which are neglected. 

If £ is of the order of 1 km and d is of the order of 1 cm, o(~) Is of magnitude, 
10’5, w hiie the second term is on the order of 2.5 x 10* 2 . Only if d is down 
to the mm level and z is down to the 0.1 km level, is the second term on the 
same order as the main term. 

A second reference system (see Figure 2-2) Is obtained by displacing 
the center of the frame to the point G. This is done when we want to consider 
the motion of the subsatellite, in the classical way, about its center of mass. 
The only torque affecting this motion in the system being analyzed here is 
the torque T due to the tension applied by the tether to the subsatellite. 

It's value is 

I = (-dk) x (-3g(0)~ a) = 3g(0)~“ k x a (3) 

Before passing to the equations of motion we define the Eulerian angle since 

different notations are used here than in the literature. 

The angle \i> is the angle measured in the plane a., j^, from a_ to N. (nodal 

line) counterclockwise (c.c.w.) with respect to an observer standing along c\ 

0 is the angle from £ to Jk measured c.c.w. from an observer oriented along ,N; 

♦ Is the angle between N and i. measured c.c.w. from an observer oriented as k^. 
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Figure 2-2. Reference Frame for Motion of Subsatellite about It's Center of Mass (G) 
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Concerning the torque given by (3), notice that w* neglect the feedback 
effect of the tether mass on the motion that we are analyzing. We also 
neglect the possible variation of the orientation of the tether from the 
vertical direction. The general computer model includes both these effects. 

2.4 Dynamical Equations 

* 

The following Is the standard procedure for writing the expression of 
kinetic energy with respect to an absolute frame 

T * j |a [e 2 + (ij>+n) 2 sin 2 e] + C [ 4 >+(f«-n)cose] 2 } (4) 

The work done by the force -JF Is given by (see Figure 2-1) 


-F • 60 ■ -j£| a_ • 6(-djc) ■ |F| da • 6,k 


■ | Fj d 6(s1nesinif<) 

« |£j d (66coses1n + 6\|is1 necosip) 

where 6 denotes a small variation In the quantity. The modules of the vector 
£ are denoted by F, below. 

The Lagrangian equations become: 


(5) 


A 

(Ae) - A(i|H-n) 2 sinecos 0 + C[^+(ij/+n)cos9] ( 4»+n) sine = Fdcosesinij< 


d_ 

dt 


^A(^+n)sin 2 e + C[$+(i{»+n)cose]cose J ■ Fdslnocos^ 
gg [c[i+(i+n)cose]j= 0 


(6-1) 

(6-2) 

(6-3) 


From the last equation we have 

♦ + (ij/+n)cos6 = r Q 

This relation means that the component of the angular momentum along the 
symmetry axis and of the angular velocity with respect to an absolute frame 
Is constant. The equations (6-1) and (6-2) become: 


( 7 ) 
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Ae - A(j>+n) 2 s1necose 4 Cr 0 (j>+n)sine ■ Fdcos0sin4 
Aislne + 2A6 (J+n)cose - CroQ * Fdcos^ 


Equations (7) and (8) give the general solution. It should be noticed, 
however, that equation (8-2) becomes singular for e * 0, and therefore, the 
numerical integration may "blow-up" at this point. 


2.5 Stationary Solution 

One may search for stationary solutions. For instance, one may search 
for solutions that satisfy the relation 

e * f 0 
o 

Equations (8-1) and (8-2) become 

-A(tj/+n) 2 sine 0 cose o + Cr 0 (^+n)s1ne o - Fdcose 0 s1m/i * 0 
Alpsi n© 0 - Fdcosi^ ■ 0 

Multiplying (10-2) by one yields 

1 Aij/ 2 s1ne o - Fdsim|> s H Q 
while equation (10-1 ) may be written 

-Aip 2 -2Aii»nsine o cose 0 - An 2 s1ne 0 cose Q 4 Cr Q («p+n)s1 ne Q - Fdcose o s1m|/ * 0 


or 


* 

Aif> 2 + Fdcose 0 sint|> - An 2 sine 0 cose 0 * 0 


and 


2Ans1ne _cose rt - Cr^slne. * 0 
oo o o 


(e-V; 

(8-Z) 


(9) 


(10-1) 

( 10 - 2 ) 


( 11 ) 

( 12 ) 

(13) 

(14) 
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For (11) to be Identic*] to (13), one should have 

7^7 ’ cos6 o’ s1n V ose o ' 2 ' 

which does not have any solution. The only solution e ■ e Q are obtained 
with ♦ ■ * 0 ■ ± 

-An 2 s1ne 0 cose 0 + Cr Q n s1ne o ± Fdcose 0 * 0 (16) 

where the signs (-) and (+) hold, respectively, when ^ ■ ± J . 

One should notice that any solution ♦ ■ j * e ■■ e Q Is coincident with 
the solution ♦ - - | , e ■ - e . Therefore, we shall consider only the case 
tl> * j . Equation (16), where we hold the sign (-) In front of the term 
Fdcose o (ij/ ■ j ) , and taking Into account (7) yields 

(C-A)n 2 s1ne 0 cose 0 + C$ 0 ns1ne 0 - "dcose o » 0 (17) 

If a o » n equation (17) reduces to 


- Fd 


“ , " D o C4> Q n 

If ci 0 and (C-A)n are comparable, we can solve equation (17). 
equation (17) In the form 


(18) 

We may write 


n 2 s1n2e 0 * 


Qs1nCe 0 -e 0 ) 


where 

Q ■ |(C 2 i 0 2 n 2 + F 2 d 2 ) 1/2 | (20-1) 

s1n6 0 * Fd/Q, cos6 0 - C* 0 n/Q (20-2) 
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Solution of equation (19) Is given graphically In Figure 2-3. The corres- 
ponding configuration* are illustrated In Figure 2-4, It Is Clear that 
when (C-A)£ « Q solution of (19) gives o Q ■ $ 6 and e Q « e 0 + (2n+l)ir. 

As an example we take: 

d ■ 0.1 m, C ■ M(0.5m) 2 , A - M(0,3m) 2 . 
t ■ 10 km, g(0) ■ 9m/sec 2 , n ■ 10" 3 , 

F ■ 3 X 9 X M X 5^5 - ^ M, 


TO” 1 . 


Then equation (17) becomes 


0.16M x 10" 6 s1ne_cose. + 0.25M 10" 4 s1ne #% 
0 0 0 

2.7M x 0.1 cose„ 


67 


( 21 ) 


which gives, In first approximation 


Tane Q • 10^, e o ■ 89?43 
e Q « 269?43 


( 22 ) 


If 4> Q * 10 ” 2 we will have 


Tane 0 * 10, e Q ■ 84® 29 
e Q » 264?29 


(23) 


Now Introduce the small angle 


A® - 6 0 - ( i Q 

where e 0 is the solution of (19), we have 


(24) 


^ n 2 (s1n2 $ 0 + 2Aecos2e 0 ) * - Qae 


(25) 


TO. 




1 


#• 


. it-'.. 
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(C-A)n 2 s1n2P 0 

Ae " ■ 2 Q + 2(C-AV cosie. 


( 26 ) 


From which we may evaluate the first order correction. 

Considering only the case f 0 ■ 10" 2 , we have for the two steady state 
configurations 

A 0 ! * - ^‘xS ' X 5l '' x^l0 3 rad “ °*^ 2 arc seconds (27.1) 

A0 2 * - 2x\ X 5 ? " xVb 3 rad “ 1,25 arc Seconds (27.2) 

We conclude by noticing that even If f 0 * O.ln, or one tenth of the mean 
motion, the correction Is negligible. 

2.6 Stability of the Equilibrium Configurations 


Since there are only two configurations of equilibrium, one presumes that 
one will be stable while the other is unstable. This may be the case, but, 
since we are confronted with the complex behavior of a three degree of freedom 
system here, other steady state solutions may exist and one must be careful 
in reaching general conclusions from this analysis in order to evaluate the 
stability of these configurations. The equations (8) are time independent 
as the Lagrangian function from which they are deduced which is: 


l » T + u = 1/2 | A [b 2 + /ij>+n) 2 sin 2 e] 
+ c[j> + (i+n)cose] 2 ^ - Fsinesini|! 


(28) 


The system has a first integral (7) which implies the invariability with 
respect to time of the component of the angular velocity vector along the 
symmetry axis, and the corresponding component of the angular momentum 

Cr Q ■ C[* + ($>+n)cose] (29) 
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Figure 2-4. Two Stationary Relative Configurations seen in the Direction of the 
Orbital Velocity Vector. 
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Moreover, the energy Integral holds 

1/2 | A e 2 + (i+n) 2 s1n 2 e + Cr 0 2 } -Fsinesinij* - E (30) 


Next we consider the linear stability of the steady state solution defined 
by the following Initial conditions: 


e(0) ■ e 0 , e(o) ■ 0, *(o) * jr. i>(0) ■ 0, 
i ■ *0 * r o • ncose o 


(31) 


where e Q is one of the two solutions of equation (16) or (17), while 4 , 
being a cyclic variable, may be choosen arbitrarily. 

We consider the equations (7) and (3) and make the following assumptions: 

6 = 0 Q + eej, Ij/ = I + eifij, iji = e$i 

♦ ■ i 0 + < 4 i» r = r Q + cr x ^ 


and substitute these in equations (7), (8), considering e a first order 
quantity. Neglecting second order quantities in e we have from (7) 

r Q + er! = 4 > 0 + 4i + ( e, i J i + 0 ) cos (e Q + eg^) (33) 

From which we get 

ri = ii + cose Q - nexsineQ (34) 

Since, however, the first integral holds also for the new initial conditions, 
we have 

ii + cos e Q - n 0 isine o = (r X ) Q (35) 

which implies that if Oj, • remain small in the pertrubed motion then 
will also remain small. Therefore, we have to study the stability with 
respect to the two variables e and ij> because the stability of $ Is 


14 . 


derived from (35). As a matter of fact, we consider only the subset of 
motions which correspond to the same value of r Q that In turn defines through 
equation (16) the value of e 0 , and ignore the $ coordinate. For other values 
of r Q different equilibrium configurations occur. For each of them the same 
assumptions and the deductions obtained here hold. 

Consider equations (8) and substitute (32) keeping r Q constant. We 
arrive at the following equations: 

eASj - |(c^ 1 +n) 2 s1n 2(e o + ee^ + Cr 0 (eii+n)s1n ( e o +e0i) 

= Fdcos (Vee^slnd- + eih) 

(36) 

ti • • • ■ f 

cA^iSln(© 0 + ce^ + 2Aee 1 (etf/ 1 +n)cos(e o + ce^ - 
= Fdcos (^ + e^i ) 

Neglecting second order quantities and taking into account that e Q is a 
solution of equation (16) we write the following linear system with constant 
coefficients 

Aei - (2Ansin0 o cose o - Cr o sin 0 o ) ^ 

+ [Cr o ncos0 o + Fdsin0 o - An 2 cos2e 0 ] ©j = 0 
A^isin0 Q + (2Anco$0 o - Cr 0 ) ©i + Fdi^ = 0 

The linear system (37) has the following characteristic equations 

(A 2 Asin0 Q + Fd)(x 2 A + Cr o ncose o + Fdsin0 Q - An 2 cos2e Q ) 

+ (2Ancos0 Q -Cr Q ) 2 sin e$x 2 

that is, the quadratic equation 


(37) 


(38) 


15 . 


ORIGINAL PAGE S3 
OF POOR QUALITY 



ORIGINAL PAGE IS 
OF POOR QUALITY 

w 2 A 2 sine o + p | AFd(l+s1n 2 e 0 ) - 3ACnr 0 cose o s1ne 0 

+ 3A 2 [n 2 cos 2 e 0 s1ne o + A 2 n 2 sin 3 e t> ]+ C 2 r Q 2 s1 ne 0 |- (39) 

+ dF (Cr„ncose_ + Fdslne - An 2 eos2e„ 

'o o o o 

where we substitute w for \ 2 . Since the system is a Lagrangian conservative 
system the characteristic exponents which are solutions of equation (39) 
either are real or imaginary; they cannot be complex. Therefore equation 

(39) certainly has real solutions and vi 2 - Stability occurs when both 
solutions are negative which means that the three coefficients of equation 

(40) should have the same sign. 

A general discussion of the signs of the coefficients are very complex. 
Some limiting cases may be treated easily. In particular, let us consider 
equation (16) choosing tf*e negative sign in front of the last term, which 
we will rewrite here, 

An 2 sine 0 cose 0 - Cr 0 nsine Q + Fdcose Q * 0 (40) 

and assume that 

An 2 << Cr 0 n» An 2 << Fd (41) 

From (40) we have 

Cr 0 nsine 0 - Fdcose Q = 0 (42) 

Similarly, equation (39) becomes 

y 2 A 2 n 2 sin 2 e 0 + usine o |F 2 d 2 cos 2 e Q + An 2 Fd (l+sin 2 e Q - 3cos 2 e Q )| 

+ FW-0 (43 ’ 
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From (43) we see that for e Q * e 0 ' the three coefficients are all 
positive and the solution Is stable, while for e o * e 0 "» the second co- 
efficient is negative and the configuration Is unstable. In the same way 
we may analyze the limiting case Fd << An 2 , Fd « Cr 0 n, and the case 
Cr R n « Fd, and Cr Q n « An 2 . 

2.7 Conclusions 

We have In this part of the report considered a very simple dynamic 
model which provides a tool for checking the numerical integration program. 
For any more complex model analytical treatment will require a larger effort 
particularly if we want to take into account possible eccentricities of the 
orbit, and the regression of the node (that is, the motion of the orbital 
plane). We could have easily found some steady state configuration for a 
triaxial body but we thought we would not have learned much more by doing 
so. However, we do believe that an extension of this analytical study to 
a more sophisticated model would be necessary for a full understanding of 
the system. In fact, while the numerical integration of the general system 
as described in this report could be quite expensive in the general case 
we think that the numerical integration of equation (6) or of the correspond- 
ing equations for a more realistic model (which may take into account 
the triaxial ity of the body, the orbital eccentricity, and the regression 
of the orbital plane) will be much less demanding In term of computing time 
and will lead to a better understanding of the dynamics of the system. 
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3.0 Rotational Dynamics Software Development and Verification 

3.1 Two-Mass Tether Model (DUMBEL) 

3.1.1 Adaptation of DUMBBELL Rotational Dynamics Model to the VAX 

Annex III of the Final Report for NASA Contract NAS8-32199, "Study of 
the Dynamics of a Tethered Satellite System (SKYHOOK)," March 1978 (attached 
as an Appendix to this report) describes a version of the DUMBEL computer 
program which models the rotation of the subsatellite. The program Integrates 
the motion of two masses - the Shuttle which is considered to be a point mass 
and the subsatellite which is modelled as a rigid-body with three moments of 
inertia and an attachment point to the tether. The rotation is described by 
twelve variables using the method of direction cosines (nine quantities) 
together with the three angular velocities about the principle axes. Six 
variables give the state vector of the center of mass of the subsatellite 
and another six give the state vector of the Shuttle for a total of twenty- 
four variables. The tether is assumed to be a massless visco-elastic con- 
nection whose dynamics is not represented by any mass points or integration 
variables. 

The rotational dynamics version of the DUMBEL program was written on 
a CDC 6400 computer and has not been used since the original study under 
which it was developed. In preparation for the current study of subsatellite 
rotation, this program was converted to the VAX computer currently in use at 
SAO. Changes have been made to allow the program to be run in an interactive 
mode with the necessary input parameters being supplied through an interactive 
terminal. (Batch mode operation can also be done using a prepared file of 
input parameters.) 

In order to verify successful conversion of this program to the VAX, an 
original CDC 6400 run was duplicated on the VAX with initial conditions of 
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\l> ■ 20°, e * 0°, and $ - 0°. All initial angular velocities were zero. The 
Initial angles are converted to direction cosines for the integration of the 
dynamics and then the direction cosines are converted back to Euler angles at 
each output point. At the beginning of the integration , the output angles 
should be the same as the input angles. After debugging* inspection of the 
original and new printouts showed only small differences between the output 
angles which are assumed to be due to the better numerical precision on the 
VAX computer In double precision. 

3.1.2 Verification of Two Mass Tether Models Using Rotational Test Case 

The test case {solved analytically in 2.0, above) consists of a special 
stationary solution of the equations of motion in which the axis of rotation 
of a symmetrical subsatellite maintains a fixed orientation with respect to 
the tether. This is accomplished by having the precession rate of the sub- 
satellite equal to the orbital angular velocity. The analytic solution was 
obtained in a reference coordinate system rotating with the orbit and used 
Euler angles to describe the orientation of the subsatellite. The axis of 
rotation is contained in the plane defined by the direction of the tether 
and the normal of the plane of the orbit. The angle between the rotation 
axis and the normal to the orbit plane is fixed and may be either plus or 
minus. Assuming the Shuttle is above the subsatellite we can define positive 
angles as measured from the normal to the orbit toward the Shuttle, and 
negative angles as measured toward the earth. As the angular velocity about 
the symmetry axis increases the angle between the rotation axis and the 
normal to the orbit decreases. For high rotation speeds, the rotation axis 
Is nearly perpendicular to the wire (and parallel to the normal to the orbit) 
in this special stationary solution. For zero rotation rate the angle to 
the orbit normal is either plus or minus 90 degrees. That is, the sub- 
satellite is hanging either straight down or straight up. The up conflgura- 
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tion is obviously unstable. With the proper initial conditions in this case, 
the subsatellite precesses at the same rate as the orbital angular velocity 
so that the subsatellite maintains a fixed orientation with respect to the 
tether. The equation giving the relation between the subsatellite angle 0 
and the spin angular velocity r is 

-A n 2 sinecose + Cnrslne ® Fdcose 


where A Is the moment of inertia in the plane of symmetry of the subsatellite, 
C is its moment of inertia about the synmetry axis (C>A), n is the orbital 
angular velocity, e is the angle between the symmetry axis and the normal to 
the orbit plane, F is the wire tension, and d is the distance from the center 
of gravity to the attachment point of the tether. For a given angle 0 we 
can solve for the required spin angular velocity r to obtain 


r 


cose / Fd 
C 'nsine 


+ An) 


The geometry must be such that the precession is in the right direction to 
make the symmetry axis follow the orbital motion. The proper geometry can 
be set up using the equation 

«9» -Jte 

N = t = d x F 

where N is the rate of change of angular momentum, t is the torque, d is 
vector to the attachment point, and f is the force due to the wire tension. 

The direction of spin of the subsatellite is shown on the circle around 
the vector d (Figure 3-1). The spin angular velocity is parallel to d and 
the torque is in the -y direction. Since the Shuttle motion is in the +y 
direction, the precession is in the proper direction. 

The parameters used in the test run are Shuttle altitude 220 km, Shuttle 
mass 86.363 metric tons, subsatellite mass 100 kg, tether length 20 km, 
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e 0 ■ 20°, * 0 * 0°, ip 0 - 90°, orbital angular Velocity n - 1.178021398 x 10" 3 
radians/sec, d ■ 50 cm, tension force F * .8342029979 x !0 6 dynes, principal 
moment of Inertia C * 2 x 10 12 , and A * 1 x 10 12 In c.g.s. units. The 
moments of Inertia have been arbitrarily chosen very large In order to obtain 
a small value of r. A fast spin rate would result In very slow numerical 
Integration. For the values of the parameters, the spin rate Is r - 

» 

-.04919339283 radian$/sec. The Initial rates of the Euler angles are * n, 

• s 

$ 0 « r, and e Q « 0. The Integration was run for 300 seconds with output 
points every 10 seconds. The angle e was .3490658 radians at the beginning 
and Increased to about .3490669 radians at the end. The value of \p was 
Initially 1.5707963 radians and increased to 1.9242006 radians at 300 seconds. 
Adding nt to gives 1,9242027 radians, The solution was stable to a few 
parts per million over the duration of the test run. 

The Euler angles of the subsatellite should be constant with respect 
to a reference frame rotating with the Shuttle. The integration is done in 
terms of direction cosines instead of angles. The direction cosines are 
converted to Euler angles at each output point. In order to check the orien- 
tation of the subsatellite in the rotating reference frame, the direction 
cosine matrix can be rotated by the negative of the orbital angle and then 
converted to Euler angles. In the test case, since the orbital motion is 
in the x-y plane, the rotated matrix is 


«1 01 Yl 

1 

cos -sin 0 

a 2 h Y2 


sin cos 0 

£*3 03 Y 3 


o 

o 

1 


where n is the orbital angle. 

When this transformation is applied in the program, the angle is constant 
during the first 300 sec to an accuracy of a few parts per million (see Figure 
3-2). 
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Figure 3-2. Spinning Subsatellite Stable Configuration* Software Test Case (a) DUMBEL: 
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3.2 Multiple-Mass Tether Model (SKYHOOK) 

3.2.1 SKYHOOK Tether Model with Subsatelllte Rotational Dynamics Capability 

The SKYHOOK computer program has been modified to Include the rotation 
of the subsatellite using the method of direction cosines as defined In 
Annex III of the report "Study of the Dynamics of a Tethered Satellite System 
(SKYHOOK)," Kalaghan et al . March 1978 (see Appendix). The angular damping 
decrlbed in Annex III has not been Included In SKYHOOK. The Inclusion of 
rotational dynamics involved changing five of the existing routines in Sky- 
hook and adding six new routines. The modified routines are the main pro- 
gram SKYHOOK, and subroutines MASSPTY , SKYIN2, DIFFUN, and TENSION. The new 
routines are DIFROT, ANGOUT, CROSS, ATTACHP, GROTAT, and PROTAT. The input 
to the program has been changed to read the number of rotational equations 
(either 12 or 0), the initial values of the rotational state vector, the 
vector from the center of mass to the attachement point of the tether, and 
the three moments of inertia about the principle axes of the subsatellite. 

The format for reading the initial conditions for the rotational 
dynamics follows the general philosophy used in Skyhook of keeping the 
input as general as possible. The user is expected to generate the initial 
conditions with a preprocessor designed for the particular problem to be 
studied. At SAO the initial conditions are generated using a small program 
called DUMBEL. The version of DUMBEL which models the rotational dynamics 
reads the Euler angles and their derivatives as the initial conditions. 

The program then calculates direction cosines and angular velocities using 
the equations on page 13 of Annex III. The initial position and velocity 
for the Shuttle and subsatellite are computed in DUMBEL from the orbital 
parameters given on input. The state vector for the Shuttle consists of 
six quantities. The total state vector for the subsatellite consists of 
18 quantities. The state vectors for each vehicle are written onto a file 
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In the format for input to the SKYHOOK program. In order to have the 
system in tension equilibrium, the offset of the center of gravity from 
the attachment point must be taken into account. This is done by computing 
the offset in the inertial coordinate system and using this vector to off- 
set the center of gravity from the attachment point so that the distance ' 

from the Shuttle to the attachment point is equal to the tether length 
specified on input. 

The following paragraphs describe the changes made to each of the 
five routines that have been changed in Skyhook. 

Three changes have been made in the main program SKYHOOK. Some 
quantities have been added to the labelled common block MASSES which occurs 
in the routines SKYHOOK, MASSPTY , SKYIN2, and DIFFUN. These quantities are 
NROTEQ (the number of rotational equations), ATTACH (3) (the vector from 
the center of mass to the attachment point in the body fixed coordinates), 

XINERT(3) (the moment of inertia about each of the principle axes), ROTMAT 
(3,3) (the direction cosine matrix), ANGVEL(3) (the angular velocity about 

each of the body axes), and ANG(3) (the Euler angles in inertial coordinates). , 

I 

The other changes involve the deployment mode of the program. If rotational } 

i 

dynamics is being included, the initial conditions for the subsatellite must 
be taken from input rather than being obtained from subroutine INITAL which 
sets up initial conditions for new masses to be deployed. 

The changes in MASSPTY consist of adding the new quantities to the 
common block MASSES and processing the additional input data needed for 
modelling the rotation of the subsatellite. The number of rotational 
equations NROTEQ is read, printed and added to the number of equations to 
be integrated for the subsatellite. Initial conditions for the subsatellite 
are read from input in the deployment mode when rotational dynamics is in- 
cluded. The point of attachment and moments of inertia are read and printed. 
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The following changes have been made to SKYIN2. The new quantities have 
been added to the MASSES common block as described in the paragraph on the 
main program SKYHOOK. Two vectors have been added providing information used 
by the numerical integrator. The vector ZMAXR(12) gives the maximum value of 
the twelve rotational quantities, and PEROT (12) gives the increment to be used 
for numerical differentiation in subroutine PEDERV. These vectors are inserted 
into the arrays ZMAX and PEDINC by a new routine PROTAT. A new subroutine 
ANGOUT is called by SKYIN2 to compute and print the Euler angles from the 
direction cosine matrix. The Euler angles have also been added to the output 
on unit 7 which is used by the plotting routines. 

Subroutine DIFFUN computes the rate of change of each variable to be 
integrated by the numerical integrator. The following changes have been made. 
The new rotational quantities have been added to the common block MASSES. The 
new routine GROTAT is called to get the rotational state vector from the master 
array. The arrays ATTACH, ROTMAT, and ANGVEL have been added to the call to 
TENSION which computes the wire tension. A new routine DIFROT is called to 
return the rate of change of the rotational quantities and the results are 
placed in the array DZ by a new subroutine PROTAT. 

Subroutine TENSION has been modified to compute the tension from the 
position and velocity of the point of attachment of the wire rather than 
from the center of mass of the subsatellite. The arrays ATTACH, ROTMAT, and 
ANGVEL have been added to the calling sequence. A new routine ATTACHP is 
called to compute the position and velocity of the point of attachment of the 
wire on the subsatellite. The tension is computed from the difference in 
position and velocity between the attachment and the neighboring mass point. 

The following six paragraphs describe the new routines that have been 
added to SKYHOOK for including the rotational dynamics of the subsatellite. 
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The rate of change of each quantity In the rotational state vector Is 
computed by subroutine DIFROT (ATTACH , ROTMAT , AN6VEL , FTENS , X INERT , DD , IM , BMASS ) , 
where ATTACH Is the offset of the attachment point of the tether from the 
center of mass of the subsatelllte, ROTMAT Is the direction cosine matrix, 
ANGVEL Is the vector angular velocity with respect to the body axes, FTENS 
Is the acceleration due to the wire tension, XINERT Is the vector giving 
the moments of Inertia about the principle axes, DD Is the output vector 
giving the rate of change of each quantity, IM Is the mass number (which 
must be 2 or the subroutine returns zero), and BMASS is the subsatelllte 
mass used to compute the wire tension from the acceleration FTENS. The 
subroutine computes the force on the subsatelllte due to the wire tension 
In the body axis coordinate system. This Is then used to compute the com- 
ponents of the torque along each body axis. The rate of change of each of 
the 12 quantities in the state vector Is computed using the equations on 
page 12 of Annex III. 

The Euler angles corresponding to the direction cosine matrix are 
computed and printed by subroutine ANGOUT ( TOUT , ZOUT , I PT , NTEQ , NEDEQ , NROTEQ , 

ANG) where TOUT is the time in seconds, ZOUT is the totat state vector, 

IPT Is a vector giving the state of each mass point in the ZOUT array, NTEQ 
Is the number of temperature equations, NEDEQ is the number of electro- 
dynamic equations, NROTEQ is the number of rotational equations for the sub- 
satellite, and ANG is the vector of Euler angles. The Euler angles are 
computed using the equations on pages 13 and 14 of Annex III. 

Subroutine CR0SS(A,B,C) computes C as the cross product of A and B. 

The position and velocity of the attachment point of the tether on the 
subsatelllte is computed by subroutine ATTACH P( ATTACH, ROTMAT, ANGVEL.SV.SVP, 

IM) where ATTACH, ROTMAT, and ANGVEL have been previously defined, S V is 
the state vector of the center of mass, SVP Is the state vector of the 



attachment point, and IM is the mass number. If IM Is not 2, SVF equals 
SV. The position of the attachment point is obtained by rotating ATTACH 
using the rotation matrix ROTMAT and adding the result to the positional 
part of SV. The velocity of the attachment point Is obtained by taking the 
cross product of ANGVEL and ATTACH and adding the result to the velocity 
components of SV. The cross product Is performed after rotating ANGVEL and 
ATTACH to Inertial coordinates. 

The rotational state vector is extracted from the total state vector 
of the system by subroutine GROTAT ( A , ND IM , I PT , NTEQ , NEDEQ , NROTEQ , ROTMAT , 

ANGVEL) where A is the total state vector for all masses, NDIM is the number 
of rows In the matrix A, IPT is the vector of indices giving the starting 
location for each mass, NTEQ Is the number of temperature equations, NEDEQ 
is the number of electrodynamic equations, NROTEQ is the number of rotational 
equations for the subsatellite, ROTMAT is the matrix of direction cosines, 
and ANGVEL is the angular velocity vector. 

The rotational state vector is placed in the total state vector array 
by subroutine PR0TAT( A, NDIM, IPT, NTEQ, NEDEQ, NROTEQ, DD) where the quantities 
are the same as described in the last paragraph for subroutine GROTAT. The 
vector DD is of length 12 and Includes both the direction cosines and angular 
velocity components. 

3.2.2 Testing of SKYHOOK with Rotational Dynamics 

The new version of the SKYHOOK program with rotational dynamics has been 
tested using the special case of Section 2.0 used previously to test the two- 
mass (DUMBEL) model. The initial rotational state vector from DUMBEL has 
been used as input to the SKYHOOK program. In testing the program It is 
helpful to plot the Euler angles to see how they change with time. The DUMBEL 
program Includes a facility for generating printer page plots of various 


quantities vs. time. The Euler angles were plotted using this facility when 
the special test case was run on the program. A similar facility exists for 
plotting the quantities written on unit 7 by the SKYHOOK program. This 
utility program has been modified to Include plots of the Euler angles which 
have been added to the SKYHOOK program output. The SKYHOOK program has been 
tested by running the test case on SKYHOOK and comparing the plots from SKY- 
HOOK with those of DUMBEL. After correction of various bugs, agreement was 
obtained between the two programs. For the test case, the angle e Is constant 
with time and the angle \ l> Increases at a rate equal to the orbital angular 
velocity (compare Figure 3-3 with Figure 3-2). The angular velocity i Is 
constant. 

3.3 Behavior of a Subsatelllte with a Single Thruster 

A simple test case has been run as an example of a rotational dynamics 
simulation with an attitude thruster. The system is assumed to be Initially 
in equilibrium with the subsatellite hanging at rest at the end of the tether 
with no. angular velocity except that of the orbital motion. At t = 0 an 
attitude thruster applies a constant torque about the body axis which is 
aligned with the tether. The parameters of the run are mostly the same as 
those used in the test case described in Section 3.1.2. For the thruster 
case the initial conditions are e 0 = ip Q = 90°, j> 0 = 1.178 x 10 rad/sec (the 
orbital angular velocity), and <j> Q = <j> 0 = 0. The Euler angles e and \p are 
shown in Figure 2-1. The torque has been chosen to be 4 x 10 8 dyne-cm 
which would result in an angular velocity of .1 radians/sec in 500 seconds 
with a moment of inertia of 2 x 10 12 (cgs). A slow acceleration allows 
the orbital a'/jgle to change significantly before the spin rate has reached 
a value which would cause slow numerical Integration. Since the equations 
of motion are written as a function of the components of the torque about 
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Figure 3-3. Spinning Subsatellite Stable Configuration. Software Test Case (a) SKYHOOK: 
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the body axes* this case can be Implemented by adding a constant to the 
torque about the body z axis. The only other torque currently modelled In 
the program Is that of the wire tension, 

The thruster case described above has been run on both the DUMBEL pro- 
gram and the SKYHOOK program. The Initial rotational state vector generated 
by DUMBEL has been used as Input to the SKYHOOK program. In each program 
the thruster model was Implemented by adding a (constant torque about the 
body z axis. The Euler angles have been plotted using the printer page as 
a graph. The output from DUMBEL differs from that of SKYHOOK In that the 
Euler angles are referred to a rotating orbital coordinate system. This 
affects only the angle <l> since the orbit Is equatorial, The difference 
between the plots from the two programs was helpful for the thruster case 
since the angle is initially constant in the rotating frame and then stops 
moving with the orbit and becomes nearly constant In the Inertial frame as 
the spin Increases. 

Figure 3-4 shows the Euler angles vs. time for the simple thruster 
case as run in the SKYHOOK program. Part a) is the inclination angle e, 
part b) is the nodal angle 4>, and part c) is the spin angle <f> about the 
principle axis. Initially, the rate of change of \p is equal to the orbital 
angular velocity. In the output of the DUMBEL program the angle was 
initially constant In the rotating orbital reference frame. As the sub- 
satellite begins to spin up, the nodal rotation is arrested as can be seen 
In part b) of the plot. Once the node stops rotating, the tether force F 
is! no longer aligned with the vector r from the center of mass to the attach- 
ment point, so there Is a torque on the subsatelllte. In Figure 3-5 we see 
that the torque T Is directed toward the +z axis. Since the angular velocity 
Is In the direction of r, the torque will cause the principle axes to move 
up out of the x-y plane, thereby decreasing the angle e as seen in part a) 
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Figure 3-4 (cont.). (c) Spin Angle 4> vs. Time. 


* 

of Figure 3-4. The angle ♦ circulates wifci: Increasing angular velocity as 
seen In part c) of Figure 3-4. The angle has not been normalized to a 
continuous variable so that there Is a discontinuity every time the angle 
passes 180 degrees. 

4.0 Adaptation of SKYHOOK General Rotational Dynamics Model to Specific 
Attitude Control Systems 

The SKYHOOK program now provides a basic tool for studying the dynamics 
of the tether system including the rotation of the subsatellite. It remains 
only to develop suitable models for "real-world" attitude control systems to 

convert their Input forces and torques 1nt®> parameters directly readable by 

* 

SKYHOOK. The implementation of such models would be done by means of routines 
that compute the torque on the subsatellite as a function of the observables 
available in the program. This torque would be added to the torque from the 
wire and the result used in the equations of motion which have been incorpo- 
rated In the program for integrating the rotation of the subsatellite. 

Routines would be developed for each of the various attitude control systems 
such as thrusters, magnetic torquers, etc. These routines are straightforward 
and do not represent a major development effort. 


39 . 




APPENDIX 


ROTATIONAL DYNAMICS OF 


SUBSATELLITE VEHICLE 


rotational dynamics of subsatellite vehicle 


m-i. equations of motion of the dumbbell system with a 

rigid body at one end 

The short dumbbell program Integrates the motion of two vehicles connected by 
A massless tether. One vehicle is treated as a point mass and the other is treated 
as a rigid body attached to the tether. Two types of damping are Included: The first 
l(s proportional to the rate of change of the distance between the two vehicles, while 
the second consists of angular damping at the point of suspension of the rigid body. 
The program has facilities for computing the initial state vector of the system from 
parameters specifying the desired orbital and system characteristics. The wire 
parameters necessary to achieve the desired initial conditions and the damping con- 
stants required for critical damping of the spring oscillations and angular oscillations 
of the end mass are also calculated. 


m-1.1 Calculation of the Initial State Vector 

A single particle in a circular orbit in an inverse-square force field obeys the 
equation 

GMm _ 2 

' ■* ' * — m rod 


or 


GM = r 3 to 2 , 


where G is the gravitational constant, M is the mass of the central body, m is the 
mass of the particle, r is the radius of the orbit, and u is the angular velocity of the 
particle. 

A similar equation can be written for two particles of mass m^ and m 2 in circular 
orbits of radius r^ and r 2 connected by a massless tether. Assuming > r^, the 
particles obey the simultaneous equations 


m-i 


and 


GMm| 



r^w 


GMm 2 



+ T * m 2 r 2 « 


2 
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which can be solved to obtain the angular velocity w and the wire tension T. We can 
eliminate T by adding the equations. This gives the equation 


m. 


( m. ».»„ \ 2 

“F + ‘X) “ <m l r l + m 2 r 2> W 
r l 


r; 


or 


m l + m 2 r 2 2 

GM = — a — — — g - 00 


(mj/rp + (m 2 /r 2 ) 


Comparing this to the equation for a single particle, we see that the system orbits 
like a single particle in an orbit of radius r given by 


-3 m l r l * m 2 r 2 

( m j/rj) + (nij/r^) 


The equation for the angular velocity is 


/5m 

VF ' 

Solving for the tension, we obtain 


T = GM 


^ r 2^ r l) “ ^ r 
(rj/mjj) + (rg/mj) 


HI-2 
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The velocities of the particles are 


and 


Vj ■ Tjtt 


v 2 * r 2 w 


The initial state vector for particles in circulsr orbits starting out on the x axis with 
Inclination I is 


*1 *1 

*l p0 > 
~ 0 , 

X 2 = r 2 

* 2 *° » 
z 2 «° r 


* 1*0 » 
yj = v x cos I , 
ij *= Vj sin I , 
x 2 -° , 

y 2 “ v 2 cosI » 

z 2 “ v 2 sin I . 


If a slight eccentricity e is desired in the orbit, we can obtain an initial state 
vector from the following simple approximation. In polar coordinates, the equation 
of an ellipse is given by 


r . a(l - e ) 

1 + e co6 (6 - 0 Q ) » 


where a is the semimajor axis. Differentiating with respect to time, we have 


a(l - e ) e sin (9 - 0 Q ) 0 


r = 


[1 + e cos (0 - 0 Q )] 4 


If e ec 1, the maximum value of r is approximately 


r mov * ae6 . 
max 
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Therefore, we can obtain a small eccentricity in the orbit by setting Xj * x 2 ■ rew in 
the initial state vector. 

m-1,2 Calculation of the Wire Constanta 

The cross-section area A of a wire of diameter D is 



If the natural length of the wire is 1 0 and the elasticity is E, the spring constant k is 



When the distance between the ends of the wire is i , the tension T is 

T ■* k(i *• f q) i l > 

= 0 , 1 < i. . 


For a dumbbell system starting out with an initial separation i and tension T, the 
natural length of the wire is calculated from 

T = k(i = -i 0 ) . 


Solving for i Q , we get 

i L 

0 1 + (T/EA) • 

HI-1.3 Critical Damping for the Dumbbell System 

The dumbbell system shown in Figure 1 oscillates about its center of mass 
(CM). 



V * 


** 
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Figure 1. 


The distances of each particle from the CM are 


t. = 


m 2 f 


1 m 1 + m 2 » 


and 


m, l 


= 


2 m 1 + m 2 


If the spring constant of the entire length is k, the spring constants of each section 


are 


and 


, m. + m 9 
k = k ~ = k — — 2. , 

1 1 1 m 2 


m l + m 2 
m. 
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■=» 

t ft The equation of motion for each section with damping la 

m /l * Vi * Vi * o » 

where bj le the damping constant. For critical damping, the determinant of the 
characteristic equation is zero, which gives 

- 4mjk } = 0 



The damping term can be rewritten in terms of f to give 



Therefore, the damping force can be written as bi, where b is defined as 
/ mm 

b 5 2 Jk \ t . 

V m i + m 2 

ni-1.4 Critical Damping for Angular Oscillations 

A simple pendulum with angular damping executing small oscillations obeys the 
equation 
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l9+b6 + Tre = 0 , 

where I is the moment of Inertia, b is the angular damping constant, T is the vertical 
force on the pendulum, r is the length of the pendulum, »nd 0 is the angular rotation 
of the pendulum from equilibrium. The characteristic equation obtained by substituting 
the solution 

8 » A e at 
is 

la? + ba + Tr = 0 , 
which has the solution 


_ -b± Vb 2 - 4ITr 
a • 

To obtain critical damping, the radical is set equal to zero, giving 

b 2 - 4ITr = 0 
or 

b = 2VTfr . 

This value of b can be used efficiently to damp out certain oscillations of the rigid 
body at one end of the dumbbell system. If the rigid body is not symmetric about 
the point of attachment, the damping constant cannot provide critical damping for 
both moments of inertia simultaneously. Motions that do not change the angle between 
r and T will not be damped. It is assumed that the last section of wire is rigid in 
order to provide a torque for the damping to work against. 

Figure 2 shows a dumbbell system with a rigid body at one end. In the diagram, 
Pj is the position of the center of mass of the rigid body, ^ is the point of attachment 
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{ si of the wire, ft is the position of the point mass, M is the mass of the rigid body, 

m is the mass of the point mass, r is the vector from to $ 2 , ^ th ® wire force 
acting on p 2 , and 6 is the angle between r and F* 



The force F is given by 

F=[k(X -i 0 ) + hl] , 

where I = |p 3 - p 2 (. The derivative of I is given by 

i= ar V<? 3 - p 2 ) • <p 3 - p 2 ) = 5 < _1 [2 (p 3 - j 2 > • ^3 ■ P2 >! 

<V? 3 > • rt* 3 - e 2 > 

( 
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The equations of motion of and tig are 

f *Mp and -ft = m $g , 

while the equation of motion for the angular momentum ia 

ft = £ , 

where £ is the angular momentum and ft, the applied torque, is 

n = n w + n d , 

in which the wire torque is ft w » r xft. The damping torque ft D is parallel to ft w 
if ft^ 4 0 or parallel to^XF+4x£if ft^ * 0. 

The applied force F is 

t = T$ , 

where T is the tension in the wire and 
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where p g = p 1 + r. 


The damping torque is 


f XF 


The quantity 6 can be obtained from differentiating either 


m-9 



or 

(dn 0 * |r X f| , 


where 



A 

F* 


t 

iil • 
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Using the cos 8 expression, we have 


at (cos = at ^ • F > f 

-0 sin 0 = $ • F + r • F , 


0 = - 


A a a A 

r • F + r • F 

* 7T 

rXF 


The above formula is not good near 6 = 0* and 180*, for which we use 
gj (sin 0) =$£ | r X F | , 


0 cos 0 


d T A A a a l 1 / 2 

= Ht[< rXF ) * ^ XF >J 


• A A 1 l 


[ A A H A A"1 

2rXF • 3 t (rXF )J > 


a 1 r X F A A A A 

ITT^Tr,- (rXF + rXF) . 


r • F rXF 


For critical damping of simple pendulum motion, the damping constant b is 
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b * 2 vTTr , 


where I le the moment of inertia, T ia the tenalon, and r ia the dietance from the 
oenter of maaa to the point of attachment. The quantity r is given by 

A , a 

r«uXr , 

where is the angular velocity of the subsatellite. The vector F is given by 


F = (a p - Ap • pgr ) , 

Ap \ |ap 1/ 


where 


Ap = p 3 -p 2 , 

• ^ a • 

Sp-P a -?2 , 

^ ^ 

= Pjj “ (P t + r ) , 

a * 

*Jh «Ja 

* P 3 - (Pj +“Xr) 


The vector r is fixed with respect to the principal axes of the subsatellite. The 
orientation of the subsatallite with respect to a set of inertial space axes can be given 
by specifying the nine direction cosines that give the transformation from the space 
axes to the body axes. The components x, y, and z of a vector with respect to the 
inertial axes are thus related to the components, x*, y r , and z r with respect to the 
principal axes of the body by the matrix equation 


m-ii 
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(*\ h a 2 

xj'k i: 

The components of the angular velocity to with respect to the principal axes are 
tip w 2 , and co 3> and the principal moments of inertia are 1^, Ig, and Ig. The com- 
ponents of L with respect to the principal axes of the body are given by the Euler 
equations 

. A*! " W 2 W 3 lI 2 “ V\ 

( ^ ■ “3“l H 3 - V ) • 

V 3 S - w i w 2 ll i • y / 

The equations for the nine direction cosines are as follows: 

+ w 2 q 3 “ w 3 a 2 = 0 * 
d 2 +U) 3 a l ’T3 =0 » 

+ w i a 2 " w 2 a l ** 0 » 

Ai+tOgPg-tOgPg^ 0 * 

^g +to 3 Pi -(0^3 = ° , 

^3 +a5 l P 2" u> 2 P l = 0 » 

\ +W 2 Y 3"V2 = 0 » 

VVl'W 0 » 

V¥2’Vl = 0 * 

giving a total of 12 equations that must be integrated to get the orientation of the sub- 
satellite as a function of time by using the method of direction cosines. 




The reverse transformation is 
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If Initial conditions are given in terms of Euler angles and their first derivatives, 
we can use the following conversion equations: 

Oj = cos ip cos <j> - cos 0 sin <}> sin ip , 

*2 - -sin ip cos <j> - cos 6 sin <j> cos ip , 
a 3 = sin 0 sin <}> , 

= cos tp sin 4> + cos 0 cos 4> sin ip , 

02 = -Bin ip sin 4> + cos 0 cos $ cos ip , 

P 3 = -sin 0 cos «}» , 

Vj = sin ip sin 0 , 

\ 2 - cos ip sin 0 , 

Y 3 - cos 0 , 

Wj, = $ sin 0 sin ip + 6 cos ip , 
o»2 = ^> sin 0 cos ip - 6 sin ip , 

Wg = | cos B + ip . 

To interpret the results in terms of Euler angles, the following equations are used: 

a 3 sin 0 sin 4> _ sin 4> 

” P 3 “ -sin 0 cos cos 9 * 

-I = sin ip sin 0 _ sin ip 
Yg ~ cos i}> sin 0 cos ip 1 

v 2 _ }/ sin 2 ip sin 2 0 + cos 2 ip sin 2 0 _ [sin o| 

Yg COS 0 = COS 0 

* 

There is an unavoidable abiguity in the sign of 0, which has been chosen to be positive 
here for convenience. 


tan 4> = 
tan ip = 

tan 0 = 
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If 9 = 0*, the angles <j> and 4> cannot be determined separately. However, the sum 
of the angles can be computed from 

P, 

tan + *l>) ■ — 

°1 

c cos (p sin 4> + cos 8 cos d> sin ip 

* cos ^ cos - cos 6 sin sin ip 

_ cos 4> sin <)> + cos <|> sin H> .« m 0 . 

~ cos tfi cob § - sin <J> sin F 4 “ 


_ sin (4> + $) 

“ cos (<}> + ifl) * 

As a simple test case, let us use a solid ball of mass M and radius r. The 
moment of inertia I is 

t _ 2 i» 2 

I = ■£ M r 


m-1.5 Test Runs of the Short Dumbbell Program with a Rigid Body at One End 

As pi-.rt of the process of testing the short dumbbell program, various cases have 
been run with different initial conditions. The program plots the tension in the wire 
plus the in-plane and out-of-plane motion of the point of attachment. 

The simplest test run was that of a subsatellite with no initial angular velocity 
and no initial angular displacement; linear and angular damping were included. An 
in-plane displacement developed initially as the wire rotated at the orbital angular 
velocity. Within about 15 sec, the subsatellite acquired the orbital angular velocity 
and followed the angle of the wire. 

Next, the subsatellite was given a 20° in-plane rotation from equilibrium with 
linear and angular damping. As the angle returned to equilibrium, the tension at first 
decreased to a minimum value in about 2 or 3 sec, increased to a maximum value at 


about 9 sec, and returned to the initial value in about 20 to 25 sec. When the same 
case was rerun, leaving out the angular damping but keeping the linear damping, the 
aubaatellite executed an angular oscillation of slowly decreasing amplitude. The 
slight damping observed results from the linear damping term, since the oscilla t i on 
causes periodic changes in the length of the wire. 

Next, the subsatellite was given a 20° rotation from equilibrium in a diagonal 
direction from the orbital plane (i. e. , both in-plane and out-of-plane displacement) 
with both linear and angular damping. A plot of the in-plane versus out-of-plane 
displacement of the point of attachment shows that it returns nearly to equilibrium 
and then executes a small amplitude rotation about the equilibrium position with con- 
stant wire tension. The same case with no damping shows an up-and-down spring 
oscillation combined with a rotational oscillation in the plane of the initial rotational 
displacement. 

In the last case run, a 20° initial diagonal rotation was combined with in-plane 
and out-of-plane initial angular velocities with linear and angular damping. The 
point of attachment ended up rotating in a small circle around the equilibrium, and the 
center of mass acquired a velocity diagonal to the orbital plane. 
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ANNEX m 

ROTATIONAL DYNAMICS OF SUBSATELUTE VEHICLE 


D. A. Arnold 



m-2. SOFTWARE DOCUMENTATION 


This program integrates the motion of two bodies connected by a massless tether. 

One body (the Shuttle) is assumed to be a point mass, and the second (the subsatellite) 

/ 

is free to rotate about the point of attachment of the wire. The program does not 
normally read any data cards. Instead, Fortran is compiled each time the program 
is run. In this way, any parameters or formulas can be changed as needed. 

Generally the subroutines remain unchanged. The following list describes the variables 
that are most likely to be changed: 


TFINAL 


number of seconds for which the motion is to be 
integrated. 

DELT 

= 

interval at which the results are to be printed (sec). 

E 

= 

elasticity of the wire material (cgs). 

DIAM 

= 

diameter of the wire (cm) . 

RLO 

= 

initial distance between the ends of the wire (cm). 

RM1 

= 

mass of the subsatellite (g). 

RM2 

= 

mass of the Shuttle (g). 

RMGAL 

s 

value of a gravity anomaly covering 1 0 X 1 ° on the 
earth's surface (mgal). 

XA, YA, ZA 

s 

coordinates of the gravity anomaly on the earth's 
surface (cm). 

HEIGHT 

s 

height of the Shuttle above the earth (cm). 

ECC 

= 

orbital eccentricity. 

B 

s 

linear damping coefficient. 

STRETCH 

= 

initial stretch of the wire (cm). 

CDRAG1 

= 

drag coefficient of the subsatellite. 

CDRAG2 

= 

drag coefficient of the Shuttle. 

AREA1 

= 

cross-section area of the subsatellite. 

AREA2 

= 

cross-section area of the Shuttle. 

RO 

= 

atmospheric-density coefficient (g/cc). 

HH 

= 

atmospheric-density scale height (cm). 

HO 

s 

reference height for atmospheric density. 

RJ2 

= 

J 2 gravitational coefficient. 
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( ; } RINC 

TttETA, PHI, PSI 
DTHETA, dphi, dpsi 
RADIUS 
XI, YI, ZI 
RB 

At each output point, the program prints the following variables giving informa- 
tion on the translational part of the motion: 

= time. 

= order of the integration polynomial. 

_ coordinates of the subsatellite and Shuttle with 
respect to the center of the earth. 

= distances of the vehicles from the center of the 
earth. 

= along-track, across-track, and radial displacements 
of the point of attachment of the wire with respect to 
the Shuttle. 

♦ 

= separation of the Shuttle from the point of attachment 
of the wire. 

= wire tension. 

= separation rate of the vehicles. 

= integration step number. 

= linear damping force. 

= integration step size. 

The angular information is printed by subroutine ANGOUT. The final state vector 
is printed and punched. The program plots the tension, along-track, and across-track 
displacements as a function of time using the printer page as a graph. It also plots 
the along-track versus the across-track displacements to give a visual display of the 
swinging motions of the subsatellite. 


TOUT 

JSTAR'i 

XI, X2A 
YI, Y2A 
Zl, Z2 j 

Rl, R2 


ALONG, ) 
ACROSS, > 
RADIAL ) 

RL 

TEN 

RLP 

NSTEP 

FP.ICT 

H 


* orbital inclination (rad)* 

« initial Euler angles of the Bubsatelllte (rad). 

* initial time derivatives of the Euler angles. 

* radius of the subsatellite (cm). 

* principal moments of inertia of the subsatellite. 

i 

= unit vector directed from the center of the subsatel- 
lite toward the point of attachment of the wire. 
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C ! m-2. 1 FUNCTION DOT (A. B) 

This routine returns the quantity X • £ as the value of the function. 

Input Parameter 

X, % = the vectors whose dot product is to be computed. 


Output parameter 
DOT 


- the dot product of A and $. 



m-i8 



m-2.2 FUNCTION XMAG (A. U) 


This routine returns the magnitude of the vector it as the value of the function, 
and a unit vector parallel to "X as the vector tT. 

Input Parameter 

it ■ any vector. 

Output Parameter 

tt = the vector £ divided by its magnitude. 

XMAG = the magnitude of the vector it. 
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m-2.3 SUBROUTINE CROSS (A. B. C) 

This subroutine computes the vector £ by the equation 

Z=tx$ , 

Input Parameter 

"X, m any vectors. 

Output Parameter 

£ = the cross product of the vectors it and 
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m- 2 .4 SUBROUTINE DIFROT <T. Y. DP) 


This subroutine computes the derivatives of the nine direction cosines and three 
components of the angular velocity. The time T is not used. 

Input Parameters 

T 
Y 

Output Parameter 

DD « derivatives of the quantities in the state vector Y. 

Subroutines called are CROSS, XMAG, and DOT. 


time (sec). 

rotational part of the state vector. 
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m-2.5 SUBROUTINE DIFFUN (T. Y. DP) 


W 


This subroutine computes the vector. DD from the state vector Y. The time T 
is not used. The forces included are the central gravitational force, the tension, if 
any, in the tether and optional perturbations, due to atmospheric drag, a gravity 
anomaly, and Jg. 

Input Parameters 

T 
Y 

Output Parameter 

DD = the derivatives of the (Quantities in the state vector Y. 

This subroutine calls SETUP, and DIFROT. 


* time (sec). 

*■ positional and angular state vector 
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m-2.6 SUBROUTINE SETUP (Y. IFLAG) 

This- ^routine does a variety of intermediate calculations necessary either for 
output or for calculating forces in subroutine DIFFUN. The results are plsced in 
common. 


Input Parameters 


Y 

* positional and rotational state vector. 

IFLAG 

* 1 if Y is a one-dimensional array 

* 2 if Y is a two-dimensional array. 

Common Variables 


RB 

* the components along the body axes of the vector 
from the center of the subsatellite to the point of 
attachment of the wire. 

RK 

* the spring constant of the wire. 

RLO 

= the natural length of the wire. 

B 

= the damping coefficient for the spring oscillations. 

GM 

*= the product of the gravitational constant and the m^B 
of the earth. 

RM1 

= the mass of the rigid body. 

RM2 

= the mass of the shuttle. 

RO 

o 

= the atmospheric density at height HO (g/cm ). 

HO 

* the height at which the atmospheric density is RO. 

HH 

* the fickle height for the atmospheric density. 

Output Parameters 


Common Variables 


RS 

= the components along the space axes of the vector 
from the center of the sub satellite to the point of 
attachment of the wire. 

DX, DY, DZ 

= the components along the space axes of the vector 
from the point of attachment of the wire on the sub- 
Satellite to the Shuttle. 

W 

= the components along the space axes of the angular 
momentum. 

DVX, DVY, DVZ 

= the components along the space axes of the velocity * 
of the Shuttle minus the velocity of the point of 
attachment of the wire on the subsatellite. 
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RL 

at 

the; distance between the ends of the wire. 

TEN 

m 

the elastic force in the wire. 

RLP 

m 

the time derivative of RL. 

FRICT 

m 

the damping force due to the rate of change of RL. 

FWIRE 

s 

the tension in the wire due to elasticity and damping. 

R1SQ 

m 

the square of the distance from the center of the 
earth to the center of the subsatellite. 

R2SQ 

K 

the square of the distance from the center of the earth 
to the Shuttle. 

R1 

« 

the square root of R1SQ. 

R2 

aa 

the square root of R2SQ. 

FG1 

= 

the central gravitational force on the subsatellite 
divided by Rl. 

FG2 

s 

the central gravitational force on the Shuttle 
divided by R2. 

V1SQ 

s 

the square of the velocity of the subsatellite. 

V2SQ 

= 

the square of the velocity of the Shuttle. 

VI 

s 

the velocity of the subsatellite. 

V2 

= 

the velocity of the Shuttle. 

CONS1 

- 

the atmospheric drag on the subsatellite divided 
by VI. 

CONS2 

- 

the atmospheric drag on the Shuttle divided by V2. 


This routine calls CROSS, which computes the cross product of two vectors. 


ni-2.7 SUBROUTINE ANGOUT (Z. XI. YI. ZI. TOUT, JSTART, NSTEP) 


This subroutine prints the angular part of the subsatellite state vector. The 
output information is the time, order of the integration polynomial, step number, 

Euler angles, a modified set of Euler angles consisting of successive rotations about 
the x, y, and z axes, the components of angular momentum along the space axes, 
the components of angular momentum about the body axes, the components of angular 
velocity about the space axes, and the components of angular velocity about the body axis. 


Input Parameters 
Z 


XI, YI, ZI 

TOUT 

JSTART 

NSTEP 


« the nine direction cosines giving the orientation of 
the body axes with respect to the inertial axes, plus 
the three components of the angular velocity with 
respect to the body axes. 

«= the principle moments of inertia. 

* time. 

*= the order of the polynomial used on the current 
integration step. 

= the step number. 


No other subroutines are called by this routine. 
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HI-2.8 SUBROUTINE ROTSTAT (THETA. Pin. PS1, PTHETA. DP HI, DPSI, Y, A) 


Tills subroutine computes the nine direction cosines and three angular velocities 
from the Euler tingles and thoir derivatives. Also computod are the components of 
the vector from the center of mass of the satellite to the point of attachment of the 
wire in the inertial coordinate system. 


Input Parameters 

THETA 

PHI 

ESI 

DTI1ETA 

DPHl 

DPSI 


* 0 . 


* 4». 

* 4k 


m i|>, 

M |J), 


0, <j>, and 4 > arc the Euler angles defining the orientation of the body axes with 
respect to the inertial axes. (5, and $ are the time derivatives of the Euler angles. 

Common Variables 


RADIES 

RB 


* the distance from the center of mass of the rigid 
body to the point of attachment of the wire, 

* a unit vootor pointing from the center of mass of 

the subsutellite to the point of attachment of the wire. 
The components are taken along the body axes. 


Output Parameters 
A 


Y 

Common Variable 
RS 


■ nine direction cosines relating the body axes to 
the Inertial axes, plus the throe components of the 
angular velocity on die body axes. 

* a matrix containing the translational state vector 
plus the rotational state vootor A. 


* the vector from the center of mass of the subsatollite 
to the point of attachment of the wire, Tho components 
are along the inertial axes. 


No other subroutines are called by tills routine. 
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